
tci1 <- function(x, alpha, basic = FALSE) {
  mean <- mean(x)
  n <- length(x)
  se <- sd(x)/sqrt(n)
  ci <- se*qt(1-(alpha/2), df = n-1)
  pval <- 2*pt(abs(mean/se), df = n-1, lower.tail = FALSE)
  
  if (basic == TRUE) {
    result <- c(mean, sd(x))
  }
  
  if (basic == FALSE) {
    result <- c(mean, mean-ci, mean+ci, pval, n)
  }
  
  return(result)
}

tci2 <- function(x, y, alpha, output) {
  n1 <- length(x)
  n2 <- length(y)
  m1 <- mean(x)
  m2 <- mean(y)
  s1 <- sd(x)
  s2 <- sd(y)
  w1 <- s1^2/n1
  w2 <- s2^2/n2
  
  v <- ((w1 + w2)^2)/((w1^2)/(n1-1) + (w2^2)/(n2-1))
  se <- sqrt((s1^2/n1) + (s2^2/n2))
  ci <- se*qt(1-(alpha/2), df = v)
  
  dom <- m1-m2
  pval <- 2*pt(abs(dom/se), df = v, lower.tail = FALSE)
  
  if (output == 1) {
    result <- c(dom, se)
  }
  
  if (output == 2) {
    result <- c(dom, dom-ci, dom+ci, pval, (n1 + n2))
  }
  
  if (output == 3) {
    result <- c(m1, m2, dom, pval, (n1 + n2))
  }
  
  if (output == 4) {
    result <- c(dom, dom-ci, dom+ci, se, (n1 + n2))
  }
  return(result)
}

bootcip <- function(x, alpha, reps) {
  props <- matrix(NA, nrow = reps, ncol = length(unique(x)))
  vals <- sort(unique(x))
  
  for (i in 1:reps) {
    draw <- sample(x, size = length(x), replace = TRUE)
    
    for (j in 1:length(vals)) {
      props[i,j] <- sum(draw == vals[j])/length(x)
    }
  }
  
  mean <- table(x)/sum(table(x))
  q <- matrix(NA, ncol = length(unique(x)), nrow = 2)
  
  for (i in 1:length(unique(x))) {
    q[,i] <- quantile(props[,i], probs = c(0+(alpha/2), 1-(alpha/2)))
  }
  
  result <- rbind(mean, q)
  row.names(result) <- c("mean", "cilb", "ciub")
  
  result <- as.data.frame(t(result))
  result$value <- rownames(result)
  rownames(result) <- NULL
  result <- result[, c(4,1,2,3)]
  
  return(result)
}

robreg <- function(spec) {
  model <- spec
  model$rob.vcov <- vcovHC(model, type="HC1")
  
  rob.pvals <- as.matrix(coeftest(model, model$rob.vcov))
  model$rob.se <- rob.pvals[,2]
  model$rob.pvals <- rob.pvals[,4]
  
  print(coeftest(model, model$rob.vcov))
  return(model)
}

boot_stats <- function(data, treat_var, treat_levels, response, wts,
                       treat_labs, alpha, its) {
  
  require(dplyr)
  require(Hmisc)
  
  if (is.null(wts)) {
    s <- select(data, treat_var, response)
    s$wts <- rep(1, nrow(s))
  }
  
  else {
  s <- select(data, treat_var, response, wts)
  }
  
  names(s) <- c("treat_var", "response", "wts")
  
  print(paste("# NA in treatment vector:", sum(is.na(s$treat_var))))
  if (sum(is.na(s$treat_var)) > 0) stop("treatment vector contains NAs")
  
  print(paste("# NA responses (dropped):", sum(is.na(s$response))))

  s <- filter(s, !is.na(response))

  fill <- rep(NA, length(treat_levels))
  tab <- data.frame(treat = treat_levels, mean = fill, cilb = fill, ciub = fill,
                    n = fill, sumwts = fill, treatlab = treat_labs)

  for (i in 1:length(treat_levels)) {

    tab[tab$treat == treat_levels[i], "mean"] <-
      wtd.mean(s$response[s$treat_var == treat_levels[i]], s$wts[s$treat_var == treat_levels[i]])

    tab[tab$treat == treat_levels[i], "n"] <- length(s$response[s$treat_var == treat_levels[i]])

    tab[tab$treat == treat_levels[i], "sumwts"] <- sum(s$wts[s$treat_var == treat_levels[i]])

    ms <- rep(NA, its)

    for (j in 1:its) {
      m1 <- filter(s, treat_var == treat_levels[i])
      m2 <- m1[sample(1:nrow(m1), nrow(m1), replace = TRUE), ]
      m3 <- wtd.mean(m2$response, m2$wts)
      ms[j] <- m3
    }

    tab[tab$treat == treat_levels[i], c("cilb", "ciub")] <-
      quantile(ms, probs = c((alpha/2), 1-(alpha/2)))
  }
  
  return(tab)
}

boot_test <- function(data, treat_var, treat_levels, response, wts, alpha, its) {
  
  require(dplyr)
  require(Hmisc)
  
  if (is.null(wts)) {
    s <- select(data, treat_var, response)
    s$wts <- rep(1, nrow(s))
  }
  
  else {
    s <- select(data, treat_var, response, wts)
  }
  
  names(s) <- c("treat_var", "response", "wts")
  
  if (length(unique(s$treat_levels)) > 2) stop("more than 2 treatment groups")
  
  print(paste("# NA in treatment vector:", sum(is.na(s$treat_var))))
  if (sum(is.na(s$treat_var)) > 0) stop("treatment vector contains NAs")
  
  print(paste("# NA responses (dropped):", sum(is.na(s$response))))
  
  s <- filter(s, !is.na(response))
  
  tab <- data.frame(estimate = NA, se = NA, pvalue = NA, cilb = NA, ciub = NA,
                    n = NA, sumwts = NA)
  
  tab$estimate <- wtd.mean(s$response[s$treat_var == treat_levels[2]],
                           s$wts[s$treat_var == treat_levels[2]]) - 
                  wtd.mean(s$response[s$treat_var == treat_levels[1]],
                           s$wts[s$treat_var == treat_levels[1]])
  
  tab$n <- nrow(s)
  tab$sumwts <- sum(s$wts)

  difs <- rep(NA, its)
  ests <- rep(NA, its)
    
    m1 <- filter(s, treat_var == treat_levels[1])
    m2 <- filter(s, treat_var == treat_levels[2])
  
    for (j in 1:its) {
      
      # For confidence interval
      
      m1_1 <- m1[sample(1:nrow(m1), nrow(m1), replace = TRUE), ]
      m2_1 <- m2[sample(1:nrow(m2), nrow(m2), replace = TRUE), ]
      
      m1_2 <- wtd.mean(m1_1$response, m1_1$wts)
      m2_2 <- wtd.mean(m2_1$response, m2_1$wts)
      
      difs[j] <- m2_2 - m1_2
      
      # For p-value
      
      rs <- s[sample(1:nrow(s), nrow(s), replace = TRUE), ]
      
      rs_m1 <- wtd.mean(rs$response[1:nrow(m1)], rs$wts[1:nrow(m1)])
      rs_m2 <- wtd.mean(rs$response[(nrow(m1)+1):(nrow(m1)+nrow(m2))], 
                        rs$wts[(nrow(m1)+1):(nrow(m1)+nrow(m2))])
      
      ests[j] <- rs_m2 - rs_m1
    }
    
  tab[ ,c("cilb", "ciub")] <- quantile(difs, probs = c((alpha/2), 1-(alpha/2)))
  tab[ ,"se"] <- sd(difs)
  
  # 2-sided p-value
  
  tab[ ,"pvalue"] <- mean(abs(ests) >= abs(tab$estimate))
  
  return(tab)
}